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1.  Introduction 

Pseudospectral  methods  have  been  used  extensively  in  simulations  of  fluid 
flows  (see,  e.g.,  Orszag  and  Patera,  1983;  Marcus,  1984;  Riley  and  Metcalfe, 
1979).  The  dependent  variables  are  expanded  in  terms  of  a  complete  set  of 
orthogonal  functions  defined  in  the  computational  domain.  Either  a  Galerkin 
method  or  a  collocation  method  can  be  employed  to  solve  the  governing  partial 
differential  equations.  The  basic  numerical  analysis  of  spectral  methods  has 
been  reviewed  in  a  monograph  by  Gottlieb  and  Orszag  (1977).  The  advantage  of 
the  pseudospectral  method  over  methods  using  locally  supported  basis  functions, 
such  as  finite  difference  or  finite  element  methods,  is  the  fast  convergence 
of  the  expansions,  which  results  in  a  numerical  scheme  with  very  small  dissipa¬ 
tive  and  dispersive  errors  for  solving  partial  differential  equations.  Appli¬ 
cation  of  the  method  to  a  compressible  flow  problem  with  shock  waves  has  been 
studied  in  one-dimensional  problems  by  Gottlieb  et  al.  (1981)  and  Taylor  et  al. 
(1981).  They  demonstrated  that  a  pseudospectral  method  can  be  used  to  capture 
shock  waves  in  the  flow  field  if  an  appropriate  filtering  technique  is  used  to 
stabilize  the  computations  and  filter  out  the  Gibb's  oscillations. 

With  the  developments  to  date,  the  pseudospectral  method  seems  to  offer 
the  opportunity  for  a  scheme  superior  to  the  finite  difference  method.  Several 
problems  need  to  be  resolved,  however,  before  a  practical  code  can  be  con¬ 
structed.  Methods  for  filtering  the  Gibb's  phenomenon  associated  with  the 
shock  wave  must  be  studied,  and  the  accuracy  of  the  results  after  filtering 
must  be  critically  examined.  The  application  of  pseudospectral  methods  to 
flows  around  a  complex  geometrical  shape  such  as  an  airfoil  has  not  been  fully 
demonstrated.  In  particular,  any  geometrical  singularity,  e.g.,  the  trailing 
edge  of  an  airfoil,  may  cause  difficulty  for  a  global  expansion  method  such  as 
the  pseudospectral  method.  The  efficiency  of  pseudospectral  computations  as 
compared  to  those  using  a  finite  difference  method  must  also  be  examined.  To 
approach  a  steady  state  through  time-marching  calculations,  a  computational 
method  allowing  large  time  steps  is  preferred.  An  implicit  pseudospectral 
scheme  may  be  a  good  candidate  for  achieving  this  goal.  A  numerical  analysis 
of  an  implicit  scheme  would  give  some  insight  into  the  efficiency  of  such  a 
scheme . 

These  are  the  issues  that  are  addressed  in  this  report.  Here,  we  present 
the  results  of  our  efforts  to  construct  a  hybrid  scheme  for  computing 
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Cransonic  flows  around  an  airfoil.  This  hybrid  scheme  consists  of  a  Fourier 
expansion  in  one  direction  and  a  finite  difference  scheme  in  the  other 
direction.  The  application  of  various  types  of  filters  is  examined  through 
numerical  experimentation.  This  work  parallels  a  recent  systematic  study  by 
Hussaini  et  al.  (1985a)  on  the  effects  of  various  filters  on  the  accuracy  of 
Euler  equation  solutions  using  a  Fourier  method  for  one-dimensional  problems. 
The  efficiency  of  our  pseudospectral  computations  as  compared  to  those  using  a 
finite  volume  scheme  is  discussed.  We  also  report  the  results  of  our  attempt 
to  construct  a  numerical  scheme  using  spectral  expansions  in  both  directions. 
Finally,  we  present  a  numerical  analysis  of  a  Richardson  iteration  scheme  for 
an  implicit  scheme.  The  analysis  is  supported  by  numerical  experiments  using 
a  simple  one-dimensional  wave  equation. 
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2.  Governing  Equations  and  Basic  Approach 

The  equations  governing  inviscid,  compressible  flows  are  the  Euler  equa¬ 
tions.  For  computing  flows  around  an  object  with  an  arbitrary  geometrical 
shape,  a  body-fitted  computational  coordinate  system  is  used.  If  (x,  y) 
represents  a  Cartesian  physical  coordinate  and  (X,  Y)  represents  the  computa¬ 
tional  space,  then  these  equations  can  be  written  in  strong  conservation  form 
by  using  contravariant  velocities  as  follows: 

It  <?>  *  lr  c?)  *  w (J)  ‘ 0  (2-1) 


where 


P 

pu 

pv 

pu 

puU  +  yyP/yM^ 

puv  -  yxP/YM^ 

q  =  J 

pv 

*► 

;  f  = 

pvU  -  XyP/W* 
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pHV 

(:)•' 

-'A 

‘x/ 

Q  ;  H  =  [(y-l)P  +  E]/p  ; 

E  =  P  + 


v(y-l) 

2 


M*  (Pu2  * 


PV2) 


(2.2) 


(2.3) 


p  is  the  density  of  the  gas;  P  is  the  pressure;  (u,  v)  are  the  Cartesian 
components  of  the  gas  velocity;  (U,  V)  are  the  unsealed  contravariant  velocity 
components  associated  with  the  curvilinear  coordinates;  J  is  the  Jacobian  of 
the  transformation;  E  is  the  total  energy;  H  is  the  specific  total  enthalpy; 

is  the  freestream  Mach  number;  and  y  is  the  ratio  between  specific  heat 
capacities.  P,  p  and  the  Cartesian  velocity  components  (u,  v)  are  nondimension- 
alized  by  their  freestream  values;  E  and  H  are  nondiraensionalized  by  the  free¬ 
stream  internal  energy  The  momentum  equations  as  written  in  the  above 

form  compute  the  evolution  of  the  Cartesian  momentum  components.  Weak  solutions 
with  shock  waves  are  allowed  by  the  above  equations,  and  entropy  increases 
across  the  shock  waves.  The  boundary  conditions  for  these  equations  are  the 
nonpermeable  conditions  on  the  solid  surfaces  and  the  assumption  that  distur¬ 
bances  generated  by  the  airfoil  are  radiated  to  infinity  with  no  reflection  at 
the  boundary  of  the  computational  domain.  For  steady-state  calculations,  it  is 
hoped  that  the  solution  will  be  insensitive  to  the  initial  conditions  specified. 
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The  basic  approach  for  applying  the  pseudospectral  method  to  flows  around 
an  airfoil  is  described  here.  The  exterior  of  an  airfoil  is  mapped  to  the 
interior  of  a  circle  using  a  conformal  mapping.  Polar  coordinates  are  used  in 
the  mapped  plane.  In  the  circumferential  direction,  we  expand  all  variables  in 
Fourier  series  because  of  the  periodicity  of  the  flow.  In  the  radial  direction 
a  Chebyshev  polynomial  expansion  for  all  variables  seems  to  be  the  natural 
choice  because  of  its  extensive  use  in  other  fluid  mechanics  calculations  in  a 
nonperiodic  domain.  In  addition  to  the  Chebyshev  expansion,  we  have  studied 
the  polynomial  subtraction  method  described  by  Gottlieb  and  Orszag  (1977)  and 
Roache  (1978).  In  the  numerical  solution  of  a  hyperbolic  partial  differential 
equation,  the  spatial  discretization  of  the  differential  operators  is  accom¬ 
plished  by  spectral  expansions.  Finite  difference  time  discretization  is  used 
to  advance  the  solution.  Because  of  their  convenience  and  computational  effi¬ 
ciency,  explicit  time-stepping  schemes  are  investigated  first.  Various  methods 
for  accelerating  the  convergence  of  the  solution  to  a  steady  state  are  imple¬ 
mented.  We  also  give  an  analysis  of  the  implicit  iterative  schemes  suggested 
by  Gottlieb  and  Orszag  (1977)  by  use  of  a  simple  one-dimensional  model  problem 
to  explore  the  numerical  properties  of  an  implicit  scheme. 
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When  a  conformal  mapping  is  applied  to  map  the  exterior  of  an  airfoil  to 
the  interior  of  a  circle,  a  shock  wave,  which  must  be  perpendicular  to  the 
airfoil  surface  for  a  transonic  condition,  will  roughly  align  with  a  radial 
coordinate  line  in  the  mapped  circle  plane.  Thus,  it  is  reasonable  to  examine 
the  capabilities  of  a  pseudospectral  method  in  resolving  a  shock  wave  by  using 
a  hybrid  scheme  that  uses  a  Fourier  series  in  the  circumferential  direction 
and  a  finite  difference  scheme  in  the  radial  direction.  The  use  of  the  hybrid 
scheme  allows  the  question  of  resolving  a  shock  wave  to  be  isolated  from  pos¬ 
sible  difficulties  that  may  arise  using  Chebyshev  series  expansions.  Techni¬ 
ques  for  filtering  Gibb's  ripples  can  also  be  studied.  For  these  reasons,  we 
have  attempted  to  construct  the  hybrid  scheme. 


3.1  Spatial  and  Temporal  Discretizations 

-f 

The  fluxes  F  and  G  in  Equation  (2.1)  are  computed  in  the  physical  space. 
Since  all  grid  points  are  defined  at  equal  intervals  in  the  transformed  space, 
the  derivatives  of  the  fluxes  along  the  circumferential  direction  are  evaluated 
by  taking  finite  Fourier  transforms  of  the  fluxes,  and  inverse  transforms  are 
applied  after  the  results  are  multiplied  by  the  wave  number.  In  the  radial 
direction,  a  second-order  central  difference  scheme  is  used  to  evaluate  all 
derivatives.  The  resulting  spatially  discretized  system  is  a  system  of 
ordinary  differential  equations  in  time. 

Following  Jameson  et  al.  (1981),  we  use  a  four-stage  Runge-Kutta  scheme 
to  advance  the  flow  variables  q  at  the  grid  points.  This  time-stepping  scheme 
given  by  the  following  procedure: 


HU)  HO)  1  jU-1).  .  .  . 

q  =  q  +  R  At;£.=  l . 4 


-*(0)  "Kn)  -Kn+1)  ■K4) 

q  =  q  ;  q  =  q 


(3.1) 

(3.2) 


where  R  is  the  residue  evaluated  by  using  q 

The  procedure  requires  only  one  level  of  memory  for  all  variables  and  is 
exactly  a  fourth-order  Runge-Kutta  scheme  for  a  linear  equation.  For  a 
nonlinear  equation,  the  procedure  does  not  have  the  high-order  time  accuracy 
of  the  true  Runge-Kutta  scheme.  However,  we  have  chosen  this  scheme  mainly 
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for  the  stability  and  not  for  the  time  accuracy.  This  simplified  four-stage 
scheme  serves  our  purpose  well.  For  a  one-dimensional  simple  wave  equation, 
the  scheme  is  stable  for  a  CFL  number  of  2.8/n  for  a  Fourier  pseudospectral 
method  as  opposed  to  2.8  for  a  finite  difference  method.  For  the  hybrid 
scheme,  the  geometric  average  of  the  time  step  given  by  the  finite  difference 
scheme  in  the  radial  direction  and  that  given  by  the  spectral  scheme  in  the 
circumferential  direction  is  used.  To  avoid  the  severe  constraint  on  the  time 
step  imposed  by  the  small  grid  size  near  the  trailing  edge,  we  use  a  constant 
CFL  number  throughout  the  entire  flow  field.  Thus,  the  time  step  at  each  grid 
point  depends  on  the  local  grid  size.  The  solution  develops  in  a  warped  time 
domain.  Again,  the  solution  will  not  be  time  accurate.  Only  the  steady-state 
solution  is  meaningful. 


3. 2  Boundary  Conditions 

It  is  well  known  that,  for  a  hyperbolic  system,  one  should  take  note  of 
various  characteristic  variables  preserved  on  the  corresponding  character¬ 
istics.  Only  the  characteristic  variables  carried  by  the  outgoing  character¬ 
istics  can  be  computed  from  the  interior  solution  of  the  governing  equations. 
The  characteristic  variables  on  the  incoming  characteristics  must  be  replaced 
by  the  appropriate  boundary  conditions.  Based  on  this  principle,  the  following 
method  of  treating  the  boundary  conditions  is  given. 

There  are  four  characteristics  crossing  the  computational  boundary.  Two 

correspond  to  the  acoustic  waves  with  wave  speed  (V  +  c),  where  V  is  the 

n  —  n 

velocity  component  normal  to  the  computational  boundary.  The  other  two  are 
the  vorticity  mode  and  entropy  mode,  both  with  wave  speed  V^.  Let  (Vg,  vn)  be 
the  flow  speeds  along  and  normal  to  the  solid  boundary,  respectively.  Then, 
the  outgoing  waves  and  the  corresponding  characteristic  variables  are 

P  -  C  V  on  -c 

c  "o  o  nc 


V  on  V  =  0 
sc  n 


P  +  yttfp  C2  on  V  =  0 
c  «  c  o  n 


(3.3 


where  the  subscript  c  denotes  the  variables  computed  from  the  interior  algo¬ 
rithm.  These  characteristic  variables  are  linearized  with  respect  to  the  flow 
variables  at  the  previous  time  step. 
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By  using  the  boundary  condition  *  0  together  with  the  characteristic 
variables  from  Equation  (3.3),  the  conservations  of  characteristic  variables 
along  the  characteristic  line  give  the  following  equations  for  the  flow 
variables  on  the  boundary  points: 


P  - 


C  V 


o  o  n 


P  - 
c 


v 


o  o  nc 


P 


(P  -  P  ) 
c 


V 

s 


=  V 

sc 


V 

n 


=  0 


(3.4) 


Other  flow  variables,  such  as  pu,  pv,  E  and  H,  can  be  recovered  from  the 
results  obtained  above. 

On  the  far-field  boundary,  flows  are  only  slightly  perturbed  from  the 
freestream  condition.  The  Riemann  invariants  are  used  to  compute  the  flow 
variables  on  the  boundary.  For  the  Riemann  invariant  on  the  outgoing  charac¬ 
teristics,  extrapolation  from  the  interior  points  is  used,  while  the  free¬ 
stream  value  is  used  for  the  Riemann  invariant  on  the  incoming  characteristics. 
After  the  normal  velocity  and  the  speed  of  sound  are  computed,  the  velocity 
tangential  to  the  boundary  is  computed  by  the  freestream  value  and  the  per¬ 
turbed  value  based  on  the  far-field  solution  of  a  circulation  around  the  air¬ 
foil.  The  rest  of  the  flow  variables  are  computed  by  the  isentropic  relations. 
We  have  attempted  other  boundary  treatments  using  a  full  characteristic  set 
similar  to  that  applied  at  the  outer  boundary.  The  results  of  computations  do 
not  seem  to  change  substantially. 


3.3  Convergence  Acceleration 

The  stability  criterion  for  the  numerical  scheme  described  is  approxi¬ 
mately  2.8/fl.  For  a  steady-state  computation,  this  requires  excessive 
computing  time.  To  increase  the  time  step  of  the  computation,  the  residue¬ 
smoothing  technique  developed  by  Jameson  and  Baker  (1983)  is  applied  to  the 
scheme  as  discussed  below. 
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Let  Che  time  change  of  the  flow  variables  be  R.  For  steady-state  computa¬ 
tions,  it  is  desirable  to  drive  R  to  zero.  Let  L  be  a  nonsingular  linear 
operator  and  R  be  defined  as 

LR  =  R  or  R  =  L_1R  .  (3.5) 

i  If  R  can  be  driven  to  zero  by  advancing  in  time,  using  R  as  the  time  change  of 

the  flow  variables,  it  automatically  ensures  that  R  becomes  zero.  Now,  R  is  a 
!  residue  defined  by  a  pseudospectra  1  spatial  expansion  of  the  variables,  then 

vanishing  R  still  ensures  spectral  accuracy  of  the  scheme  even  if  L  is  not  a 

spectral  operator.  The  choice  of  L  is  to  improve  the  stability  of  the  numeri- 

i 

;  cal  scheme.  It  has  been  shown  by  Jameson  and  Baker  that  the  following 

|  operator  serves  the  purpose: 

j  L  *  (1  -  e62)  (1  -  e6^)  (3.6) 

where  6^  and  6^  are  finite  difference  operators  in  the  transformed  coordinates 
[  (X,  Y).  To  understand  the  reason  for  improved  stability  for  a  pseudospectral 

|  scheme,  a  simple  one-dimensional  wave  equation  is  considered: 

<J>t  +  c<l>x  “  0  •  (3.7) 


The  residue-smoothing  process  as  described  is  equivalent,  to  the  lowest  order, 
to  adding  an  additional  term  to  the  original  simple  wave  equation  and  convert¬ 
ing  it  to  the  following  equation: 

+  cd>  -  e(Ax)^$.  ■  0  .  (3.8) 
Tt  Tx  txx 


The  dispersion  relation  for  this  equation  can  be  given  as 

h)  m  _ c _ 

k  1  ♦  ek2(Ax) 2 


(3.9) 


where  U)  is  the  frequency  and  k  is  the  wave  number.  By  increasing  the  para¬ 
meter  £,  the  wave  speed  for  the  high- wave -number  component  is  substantially 
reduced.  This  reduction  in  wave  speed  for  the  dangerous  short  waves  contri¬ 
butes  to  the  increase  in  the  time  step  size.  Equation  (3.8)  is  the  linearized 
form  of  a  model  equation  for  long  dispersive  waves  proposed  by  Benjamin  et  al. 
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(1972).  They  pointed  out  the  numerical  advantage  of  using  this  equation  over 
the  Kortweg-de-Vries  equation.  Other  means  of  manipulating  the  dispersion 
relation  to  gain  stability  have  been  suggested  (see  Gottlieb  and  Turkel,  1980, 
for  example).  However,  these  methods  do  not  preserve  the  original  conserva¬ 
tion  laws  as  their  limit  in  steady  state. 

In  principle,  the  method  described  above  can  be  stable  for  any  desired 
time  step  if  a  large  enough  £  is  chosen.  In  practice,  the  approach  to 
steady  state  can  be  delayed  if  a  large  £  is  chosen.  We  take  E  *  1.5,  and 
the  CFL  can  be  increased  50%  from  its  original  value.  Beyond  this  value,  it 
does  not  seem  to  help  the  approach  to  steady  state. 

3.4  Filter  and  Artificial  Viscosity 

The  above  scheme  by  itself  is  unstable  due  to  the  buildup  of  high-frequency 
components  of  error.  In  the  finite  difference  scheme,  a  form  of  artificial 
viscosity  has  been  suggested  by  Jameson  et  al.  (1981)  for  incorporation  into  a 
central  difference  scheme.  Sakell  (1984)  has  chosen  second-order  and  fourth- 
order  viscosity  in  his  study  of  pseudospectral  schemes.  It  is  well  known  that 
the  use  of  artificial  viscosity  will  broaden  the  shock  thickness.  The  objec¬ 
tive  of  using  a  pseudospectral  method  to  better  resolve  the  shock  wave  is  thus 
defeated.  Furthermore,  we  have  experimented  with  the  artificial  viscosity 
using  the  form  given  by  Jameson  et  al.  (1981).  We  have  found  that  the  amount 
suggested  by  Jameson  et  al.  is  not  sufficient  to  control  the  Gibb's  phenomena, 
and  the  amount  must  be  increased  to  obtain  a  stable  calculation.  This  increase 
further  broadens  the  shock  wave  to  an  unacceptable  level.  For  this  reason,  we 
have  decided  to  use  the  artificial  viscosity  in  Jameson's  form  only  in  the 
radial  direction,  along  which  a  finite  difference  discretization  has  been 
used.  In  the  circumferential  direction,  another  kind  of  filter  must  be  used. 

We  have  conducted  a  systematic  experiment  with  various  types  of  filters 
for  pseudospectral  transonic  calculations.  Our  experience  is  summarized  here. 

High  Wave-Number  Cutoff 

After  advancing  a  time  step,  the  flow  variables  are  expanded  into  Fourier 
series.  The  coefficients  of  the  highest  10%  wave  numbers  are  set  to  zero.  It 
is  found  that  the  time  marching  cannot  be  stablized.  The  intermediate  solu¬ 
tions  also  show  that  the  shock  resolution  has  deteriorated.  The  method  is 
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equivalent  to  adding  artificial  viscosity  terms  of  even  order  to  the  conserva¬ 
tion  equations.  Other  filters  in  wave-number  space  have  similar  properties 
and  do  not  stabilize  the  solution.  Similar  conclusions  were  obtained  be 
Hussaini  et  al.  (1985).  A  point  worth  mentioning  is  that  the  solution  in  the 
present  scheme  is  advanced  in  the  physical  space  as  opposed  to  advancing  the 
Fourier  coefficients  in  time.  Filtering  in  the  spectral  space  requires 
additional  finite  Fourier  transforms  of  all  flow  variables.  The  operation 
count  of  the  computation  is  substantially  increased. 

Derivative  Filtering 

As  will  be  discussed  later  (Section  3.5),  the  residue  of  the  pseudospectra  1 
scheme  does  not  decrease  monotonically  during  time  marching  for  steady-state 
computations.  One  of  the  reasons  is  the  Gibb's  error  in  evaluating  spectrally 
the  derivatives  of  the  fluxes.  Consider  a  step  function  in  physical  space. 

The  derivatives  are  zero  everywhere.  By  using  finite  Fourier  transforms  to 
evaluate  the  derivatives,  a  large  error  is  comnitted  near  the  discontinuity. 

We  have  attempted  to  filter  the  derivatives  of  the  fluxes  in  the  spectral 
space  when  forming  residue  in  the  hope  of  obtaining  a  better  estimate.  The 
numerical  experiments  show  that  this  type  of  filtering  cannot  stabilize  the 
computations. 


One-sided  Schumann  Filter 

Another  form  of  filtering  is  to  average  the  flow  variables  around  the  grid 
point.  In  one  dimension,  this  takes  the  following  form: 


q.  -  0.25  (q^  ♦  2q£  +  q^)  .  (3.10) 

To  avoid  smearing  shock  waves,  the  sonic  points  at  the  shock  are  located  on 
each  coordinate  line  around  the  airfoil.  The  following  one-sided  averaging  is 
applied  on  each  side  of  the  shock  wave. 

q£  -  0.5  (q.  ♦  q.+1)  (3.11) 

This  filter  has  been  suggested  by  Gottlieb  et  al.  (1981).  If  applied  at  each 
time  step,  it  is  equivalent  to  a  first-order-accurate  second-order  artificial 
viscosity.  We  have  applied  this  filtering  once  every  L  time  steps,  just  often 
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enough  to  damp  the  instability.  It  is  believed  that  this  application  of  the 
filter  reduces  the  error.  This  filter  has  been  found  to  be  the  most  effective 
in  stabilizing  the  calculations.  The  amount  of  artificial  viscosity  is 
roughly  of  the  order  of 


(Ax)2 

4(LAt) 


(3.12) 


Assuming  that  LAt  is  of  0(1)  as  in  the  present  investigation,  the  artificial 
viscosity  is  of  the  order  of  (Ax)  in  contrast  to  (&)  for  the  form  suggested 
by  Jameson  et  al.  (1981).  As  has  been  discussed,  the  form  of  artificial  vis¬ 
cosity  suggested  by  Jameson  et  al.  is  not  sufficient  to  stabilize  the  pseudo- 
spectral  computation.  It  seems  that  a  lower  order  filter  such  as  the  present 
one  is  required.  The  advantage  of  the  present  filter  is  that  the  sharpness  of 
the  shock  wave  can  be  maintained. 


3.5  Computed  Results 

Test  calculations  were  performed  for  the  flow  around  a  Karman-Tref ftz 
airfoil.  An  analytic  transformation  that  maps  the  interior  of  a  circle  in  the 
computational  space  £  to  the  exterior  of  an  airfoil  in  the  physical  space  Z 
is  given  as  follows: 

Z  ,  (1+U)K  ♦  (l-Jk)* 

Kl  (1 +££)*  -  (l-iO'C 

with 

i  -  (i  -  n‘)l/2  -  S0 

where  *  (£Q,  n  )  is  the  coordinate  of  the  center  of  the  circle  in  the 
C  plane.  For  the  present  example,  *  (-0.1,  0).  At  the  trailing  edge,  this 
transformation  is  singular.  In  the  present  investigation,  the  trailing  edge 
is  placed  at  the  half  grid  between  grid  points  to  avoid  the  singularity.  The 
coordinates  of  the  grid  points  in  the  physical  space  are  computed  from  the 
transformation,  and  the  transformation  matrices  are  computed  at  each  grid 
point  using  exactly  the  same  algorithm  for  the  evaluation  of  flux  derivatives. 
These  matrices  are  stored  to  save  computational  time.  Figure  1  shows  the 
computational  grid  around  the  airfoil. 


(3.13) 


(3.14) 
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The  computations  start  from  an  incompressible  velocity  field  determined  by 
the  conformal  mapping.  The  density  and  pressure  fields  are  given  by  the  isen- 
tropic  relations.  The  lift  and  drag  coefficients  are  computed  by  integrating 
the  pressure  around  the  airfoil.  To  be  consistent  with  the  pseudospectral 
numerical  scheme,  the  integrals  have  been  performed  by  using  a  Fourier  series 
expansion  of  the  surface  pressure  as  follows: 


Let 


-(^)P  n  da  -  ~(^)|p  y: 


XI  t  “  F  *y 
**Y«0  |Y 


..51 


dX 


(3.15) 


*P  yx|  "  fx 

lY*0 

P  xY|  -  f 

i  Y“0  y 


(3.16) 


Then  the  force  F  ■  (F  ,  F  )  is  simply  the  coefficients  a  and  a  of  the 

x  y  ox  oy 

Fourier  expansion  for  f  and  f  . 

r  x  y 

A  transonic  case  with  a  frees tream  Mach  number  of  0.63  and  an  angle  of 
attack  of  2  degrees  is  computed.  The  computations  are  done  on  both  a  64x16 
grid  and  a  32x16  grid.  Figure  2  shows  the  initial  distribution  of  the  pressure 
coefficient  on  the  denser  grid,  and  Figure  3  shows  the  solution  after  2400  time 
steps.  The  pressure  distribution  essentially  remains  constant  after  1800  time 
steps.  Figure  4  shows  the  results  of  a  computation  using  the  finite  volume 
code  FL053  with  the  same  grid  density.  Notice  that  a  weak  shock  wave  appears 
on  the  pseudospectral  calculation.  The  same  case  is  computed  on  a  32x16  grid, 
and  the  results  are  shown  in  Figure  5.  The  results  indicate  that,  although 
the  32x16  grid  gives  the  gross  features  of  the  pressure  distribution,  the 
solution  is  not  accurate  enough.  The  drag  count  for  the  spectral  calculations 
seems  to  be  too  large  for  this  low  supercritical  case,  which  may  be  due  to  the 
error  introduced  by  the  filtering. 

A  supercritical  case  with  a  freestream  Mach  number  of  0.7  and  an  angle  of 
attack  of  2  degrees  has  also  been  computed  on  the  64x16  grid.  Figures  6 
through  8  show  the  pressure  distributions  after  0,  800,  and  2400  time  steps, 
respectively.  The  solution  seems  to  reach  a  steady  state  in  1600  time  steps. 
This  is  approximately  equivalent  to  the  500  time  steps  required  in  the  finite 
volume  computation  if  we  take  into  account  the  ratio  of  IT  between  the 
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time  steps  for  the  two  methods  to  satisfy  the  stability  requirements.  The 
shock  wave  is  resolved  in  one  grid.  Figure  9  shows  the  isomach  lines  in  th“ 
flow  field.  Figures  10  and  11  show  the  results  from  a  finite  volume  calcula¬ 
tion  on  the  same  grid.  The  shock  resolution  is  not  as  good  as  in  the  spectral 
calculation.  We  have  performed  a  finite  difference  calculation  with  a  central 
difference  scheme  in  the  circumferential  direction.  By  using  the  same  filter 
as  in  the  pseudospectral  calculation,  we  intend  to  determine  whether  the  shock 
resolution  is  the  result  of  the  particular  filter  or  the  result  of  pseudospec¬ 
tral  discretization.  The  finite  difference  calculation  diverges,  which  seems 
to  indicate  the  merit  of  the  spectral  scheme.  Figure  12  shows  the  results  of 
a  32x16  computation  using  the  hybrid  scheme.  Again,  it  gives  the  correct 
gross  features  of  the  flow  field,  while  the  finite  volume  scheme  on  the  same 
grid  gives  unacceptable  results  with  a  smeared  shock  wave. 

The  results  of  the  computations  can  be  summarized  as  follows: 


(1)  The  shock  wave  can  be  resolved  within  one  grid  using  the  hybrid  scheme, 
while  it  takes  three  to  four  grids  for  the  finite  volume  code.  The 
best  upwind  finite  difference  codes  (Anderson  et  al.,  1985;  Yee, 
private  communication,  1985)  can  resolve  a  shock  wave  within  two  grids. 


(2)  The  accuracy  in  the  smooth  flow  region  seems  to  be  deteriorated  by  the 
filter  being  used.  The  expansion  on  the  lower  surface  is  probably  not 
accurate  enough,  and  the  drag  count  is  too  high. 


(3)  The  residue  in  computations  does  not  decrease  with  time.  The  residue 
decreases  first  and  then  starts  to  diverge.  At  this  point,  the 
Schumann  filter  is  applied,  and  the  residue  resumes  its  decreasing 
trend.  The  computation  is  a  slowly  divergent  one  stabilized  by 
periodic  filtering.  Hence,  it  is  very  difficult  to  decide  the 
steady-state  solution.  In  fact,  for  supercritical  cases,  the 
variation  of  the  number  of  supersonic  points  is  used  as  an  indicator 
for  steady  state. 


(4)  The  computations  with  the  pseudospectral  scheme  require  more  time 
steps  to  reach  a  steady  state  due  to  the  smaller  time  steps  required 
by  the  CFL  restriction.  To  compute  the  supercritical  lifting  case 
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presented  in  the  present  report,  1800  time  steps  are  required  as 
opposed  to  500  time  steps  for  the  finite  volume  scheme  with  the  same 
number  of  grid  points.  The  ratio  is  roughly  It  as  indicated  by  the 
stability  analysis.  The  additional  number  of  time  steps  can  be 
justified  only  if  the  pseudospectral  calculations  can  reduce  the 
number  of  grid  points  by  a  factor  of  4.  Another  possibility,  which 
was  not  investigated  here,  is  to  use  a  multigrid  method  in  conjunction 
with  the  pseudospectral  method  to  reduce  the  number  of  work  cycles. 

These  examples  show  that  the  pseudospectral  method  may  not  be  competitive 
with  the  finite  difference  methods  for  steady-state  problems  with  shock  waves. 
The  ability  of  the  pseudospectral  calculations  to  resolve  a  shock  wave  in  one 
grid  does  not  justify  the  excessive  computing  time  required.  Spectral  accuracy 
is  destroyed  by  the  existence  of  discontinuity.  Since  the  Gibb's  error  is  a 
systematic  error,  one  would  hope  that  the  spectral  accuracy  can  be  restored  by 
a  proper  filtering  technique.  The  filters  examined  in  the  present  work  fail 
to  accomplish  this  purpose. 
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Development  of  a  Fully  Spectral  Method 


In  this  section,  we  describe  the  results  of  an  effort  to  extend  the  method 
to  include  a  spectral  expansion  in  the  radial  direction  in  the  computational 
space.  If  successful,  the  scheme  will  be  a  fully  spectral  scheme  with  high 
accuracy  in  both  directions. 

Chebyshev  polynomial  expansions  have  been  extensively  used  in  the  simula¬ 
tions  of  incompressible  viscous  flows  (see,  e.g.,  Orszag  and  Patera,  1983; 
Marcus,  1984).  Application  of  the  Chebyshev  method  to  compressible  flows  using 
the  full  potential  equation  has  been  investigated  by  Streett  (1983).  Hussaini 
et  al.  (1985b)  used  this  method  in  conjunction  with  a  shock-fitting  method  to 
study  the  interactions  between  a  shock  wave  and  a  hot  spot.  They  have  shown 
that  the  Chebyshev  method  is  highly  accurate  in  that  case.  The  application  of 
the  method  to  a  steady-state  compressible  flow  using  a  shock-capturing,  Euler 
equation,  time-marching  technique  has  been  investigated  by  Gottlieb  et  al. 
(1984).  It  seems  to  be  a  natural  choice  for  developing  a  fully  spectral  scheme 

A  polynomial  subtraction  Fourier  method  has  been  suggested  by  Gottlieb  and 
Orszag  (1977).  Roache  (1978)  studied  the  accuracy  of  the  method  in  conjunction 
with  the  degree  of  the  polynomial  used.  Orszag  and  Patera  (1983)  used  the 
method  for  simulations  of  incompressible  viscous  flows.  We  shall  investigate 
the  possible  application  of  the  method  to  the  present  case. 

4.1  Fourier  Chebyshev  Method 

In  the  unit  circle  computational  plane,  the  dependent  variables  are 
expanded  into  a  finite  Fourier  series  in  the  circumferential  direction  and  a 
Chebyshev  series  in  the  radial  direction.  In  the  radial  direction,  the  compu¬ 
tational  domain  (e,  1)  is  divided  into  (N+l)  unequal  Chebyshev  collocation 
points.  In  the  present  work,  the  collocation  method  is  used. 

Derivatives  of  Fluxes 

The  fluxes  are  computed  at  the  grid  points.  The  radial  derivatives  of  the 
fluxes  are  computed  by  a  finite  Chebyshev  expansion  of  these  data.  Due  to  the 
nature  of  the  transformation,  the  fluxes  are  singular  at  the  center  of  the 
circle.  The  singularity  behaves  as  follows: 

(F,  G)  “  K/Y  (4.1 
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where  K  is  the  strength  of  the  singularity.  Although  the  far-field  boundary  is 
chosen  at  a  finite  distance  Y  =  0.02,  substantial  errors  are  committed  in  the 
spectral  evaluations  of  the  derivatives  of  fluxes  due  to  the  inability  of  the 
finite  series  to  sufficiently  resolve  the  singularity  if  17  or  fewer  grid 
points  are  used.  As  a  result,  ripples  occur  in  the  evaluated  derivatives. 

These  ripples  are  not  confined  to  the  far  field  but  also  propagate  to  the  near 
field.  To  reduce  these  ripples,  we  have  subtracted  a  singular  function  from 
the  radial  flux  vectors  as  follows: 

£  -  G  -  it/Y  (4.2) 

where  K  is  determined  by  the  fluxes  at  the  first  two  points  in  the  far  field. 

The  reduced  fluxes  G'  are  expanded  in  Chebyshev  series.  The  spatial  deriva¬ 
tives  of  G  are  obtained  by  the  combination  of  the  derivatives  of  G'  and  the 
singular  term.  This  method  has  proven  to  eliminate  the  ripples  in  the 
evaluation  of  derivatives  by  Chebyshev  series.  However,  the  process  increases 
the  computational  overhead. 


Stability  and  Convergence  to  Steady  State 
The  following  formula  defines  the  Chebyshev  collocation  points  for  an 
N-term  series: 


Y. 

J 


(1  -  eo) 


1  -  cos 


*(j-l) 


N 


+e;j*l,  ...,  N+i  .  (4.3) 


,-2, 


These  points  are  close  together  near  the  boundary  with  spacing  of  0(N  ). 

The  stability  requirement  is 

r 

1 


At  -  0 


NZ(U+C) 


(4.4) 


where  the  characteristic  velocity  is  (U+C).  The  computational  time  required 
to  reach  a  steady  state  is  increased  by  a  factor  of  N  as  compared  to  that 
using  an  equally  spaced  Fourier  method.  This  is  in  contrast  to  the  incompres¬ 
sible  viscous  simulations,  where  the  characteristic  velocity  approaches  zero 
on  the  boundary.  The  time  steps  can  maintain  a  reasonable  value  even  with  a 
highly  stretched  grid,  like  the  Chebyshev  collocation  grid  in  that  case. 

The  residue-averaging  technique  and  the  use  of  local  time  steps,  as 
described  in  the  previous  section  for  the  hybrid  scheme,  fail  to  improve  the 
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convergence.  To  understand  the  performance  of  the  Chebyshev  method  in  a  com¬ 
pressible  flow,  particularly  with  the  boundary  conditions,  we  have  tested  a 
wave  propagation  problem  in  one  dimension.  The  Euler  equations  are  solved  in 
one  dimension  with  a  smoothed  pressure  discontinuity,  isothermal  and  no  flow 
as  initial  conditions.  The  pressure  difference  across  the  wave  is  2.  We  are 
interested  in  the  time  step  size  and  the  effects  of  boundary  conditions  on  the 
stability  of  the  computations.  Figure  13  shows  the  results  of  computation 
using  33  grid  points.  In  order  to  maintain  stability,  the  time  step  is  so 
small  that  it  takes  4800  time  steps  for  the  pressure  wave  to  propagate  across 
half  the  computational  domain.  The  boundary  treatment  does  produce  shock 
reflection  on  the  solid  boundary.  The  main  conclusion  obtained  here  is  that 
unless  an  acceleration  scheme  can  be  found  to  drive  effectively  the  solution 
to  steady  state,  the  Chebyshev  explicit  method  is  not  suitable  for  steady- 
state  calculations.  Since  the  use  of  a  local  time  step  and  residue  averaging 
fails  to  stabilize  the  calculations,  an  implicit  method  may  be  the  solution  to 
the  problem.  This  is  the  reason  we  have  decided  to  investigate  an  implicit 
iterative  scheme,  as  described  in  Section  5. 


4.2  Polynomial  Subtraction  Method 

To  circumvent  the  difficulty  of  the  Chebyshev  method  as  discussed  above, 
we  have  experimented  with  a  polynomial  subtraction  Fourier  method. 

Let  G  be  the  flux  vector.  The  computational  space  ( e,  1)  is  divided  into 
N  equally  spaced  grids.  Consider  a  reduced  flux  vector  G'  defined  as 


G'  =  G  -  K/Y  -  P  (Y)  (4.5) 

m 

where  K  is  the  strength  of  the  singularity  at  Y  =  0  as  defined  in  the  previous 
section  and  Pm  is  a  polynomial  of  mth  degree.  One  can  construct  a  polynomial 
Pffl  in  such  a  way  that  G'  is  periodic  at  the  boundaries  up  to  the  (m-l)th 
derivatives  at  its  boundaries.  The  derivatives  of  the  flux  vector  G  can  be 
reconstructed  by  using  the  Fourier  spectral  method  on  G'  and  analytical 
expressions  for  the  other  two  terms. 
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Accuracy  of  the  Method 

Since  the  method  only  eliminates  the  discontinuity  of  G  up  to  the  (m-l)th 
order,  the  resulting  Fourier  series  converges  as  For  example,  if  a 

third-order  polynomial  is  used,  the  first  discontinuity  at  the  boundaries  is 
the  third  derivative.  The  resulting  Fourier  method  is  fourth-order  accurate. 
This  estimation  is  based  on  the  assumption  that  the  'jumps'  of  the  derivatives 
of  G  at  the  boundaries  are  known  exactly.  In  the  present  problem,  these 
'jumps'  must  be  evaluated  by  approximate  means.  The  error  committed  in  the 
evaluation  of  the  jumps  must  be  consistent  with  the  order  of  the  polynomial 
being  used. 


Construction  of  Polynomial  Pm 

In  order  to  construct  P^,  the  jumps  in  the  derivatives  of  the  flux  vector 
G  up  to  the  (m-l)th  order  at  both  ends  of  the  finite  domain  are  required. 

Since  G  is  known  only  at  the  discrete  points,  the  evaluation  of  these  deriva¬ 
tives  can  only  be  achieved  by  using  a  finite  difference  method.  The  truncation 
error  in  this  process  will  also  affect  the  accuracy  of  the  method.  For  example, 
when  a  second-order  one-sided  differencing  is  used  to  evaluate  the  end  point 
derivatives,  the  jump  in  the  first-order  derivative  of  the  reduced  flux  G’  at 
the  boundary  is  of  0(Ax  ).  Since  the  0(1)  jump  in  the  first-order  derivative 
results  in  0(N  )  convergence  in  the  Fourier  series,  the  resulting  polynomial 

subtraction  method  is  of  0(N  ^).  Hence,  it  is  consistent  to  use  a  second-order- 
accurate  evaluation  of  the  jump  in  conjunction  with  the  use  of  a  third-order 
polynomial. 


Computational  Experiments 

We  have  experimented  with  the  method  using  a  third-order  polynomial.  A 
computation  with  the  same  supercritical  case  computed  earlier  using  the  hybrid 
scheme  is  carried  out.  The  time  step  for  the  'full  spectral'  calculations  is 
half  that  for  the  hybrid  scheme.  The  computations  produce  similar  results  to 
the  hybrid  scheme,  with  computational  time  more  than  double  that  of  the  hybrid 
scheme.  We  have  decided  that  the  method  is  not  competitive  with  a  hybrid 
scheme  with  twice  the  number  of  grid  points  in  the  radial  direction. 
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5.  Implicit  Pseudospectral  Method 

There  are  two  reasons  why  we  would  like  to  examine  an  implicit  method. 

First,  as  discussed  in  the  previous  sections,  there  are  several  problems  in 
using  explicit  schemes  for  pseudospectral  Euler  equation  calculations.  One  of 
the  most  difficult  problems  with  these  schemes  is  the  6mall  time  steps  required 
by  the  stability  criteria.  The  Chebyshev  method,  in  particular,  requires 
extremely  small  time  steps.  An  implicit  method  may  free  us  from  this  con¬ 
straint.  The  other  reason  is  to  construct  an  implicit  method  for  time-accurate 
unsteady  calculations.  Consider  the  flow  around  an  airfoil.  The  grid  size 
near  the  leading  and  trailing  edges  is  extremely  small.  For  an  unsteady 
calculation,  one  cannot  use  the  local  time  step  as  described  in  Section  3.  An 
implicit  method  is  essential  in  this  case  to  avoid  excessive  computing  time. 

Table  1  summarizes  the  performance  of  various  pseudospectral  schemes  for 
the  wave  equation 

u  +  u  =0  (5.1) 

t  x 

4 

It  can  be  seen  that  an  operation  count  of  0(N  )  is  required  for  all  implicit 
methods  if  a  direct  inversion  method  16  used.  Since  the  pseudospectral  method 
gives  a  full  matrix  operation  on  the  variables  at  grid  points,  a  sparse  matrix 
technique  for  inverting  the  operator  is  not  applicable. 

We  shall  examine  an  iterative  implicit  method  that  uses  a  finite 
difference  operator  as  an  approximate  operator  for  the  spectral  operator  for 
the  iteration.  If  there  exists  a  stable  and  convergent  iterative  scheme,  the 
work  required  for  such  a  scheme  will  be 


S 


K 


Is 


W  =  MN  log  N 


(5.2) 


where  N  is  the  number  of  modes  and  M  is  the  number  of  iterations  per  time  step. 
Here,  we  have  assumed  that  NAt  “  0(1)  for  time  accuracy,  and  for  each  itera¬ 
tion  the  operation  count  is  roughly  that  required  by  an  explicit  scheme,  i.e., 
the  overhead  on  inverting  an  approximate  operator  is  neglected.  Therefore, 
the  Chebyshev  implicit  iterative  method  will  have  better  performance  than  the 
explicit  method  if  M  is  much  less  than  0(N).  For  a  steady-state  calculation 
where  it  is  possible  to  have  NAt  >>  1,  the  method  can  also  improve  the 
performance  of  a  Fourier  scheme. 
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Table  1.  Estimated  Performance  of  Various  Pseudospectral  Schemes 


Fourier 

Leap  Frog 

Fourier 

Crank-Nicolson 

Fourier 

Forward  Euler 

Fourier 

Backward  Euler 

Fourier 

2-step  Runge-Kutta 

Fourier 

3-step  Runge-Kutta 

Fourier 

4-step  Runge-Kutta 

Chebyshev 

Leap  Frog 

Chebyshev 

Adams-Bash forth 

Chebyshev 

Crank-Nicolson 

Chebyshev 

Backward  Euler 

A 

Tt 

2iN  log  N 

Stable 

N4/3 

N3 

C/N 

—  log  N 

Stable 

N4/3 

v  3 
u 

2 tt/T  N2  log  N 

CM 

2  TT  yfz  N2  log  N 

Unstable 

— 

4/N 

n3  1  N 

—  log  N 

Stable 

N4/  3 

Stable 

N4/3 
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We  shall  examine  a  Richardson  iteration  scheme  in  this  section.  Two 
problems  regarding  the  iterative  scheme  are  addressed  here:  What  is  the  best 
way  to  construct  a  scheme  that  will  converge  quickly  at  each  time  step?  Is 
the  time-stepping  scheme  stable  if  only  a  finite  number  of  iterations  at  each 
time  step  are  carried  out?  These  problems  are  addressed  through  the  numerical 
analysis  of  a  simple  wave  equation  and  numerical  experiments  to  substantiate 
the  theory. 


5.1  Richardson  Iteration 

Consider  the  simple  wave  equation 


3u  3u 
3t  ax 


0 


(5. 


In  finite  dimensional  space,  we  approximate  the  above  equation  by  the 
following  backward  Euler  pseudospectral  scheme: 

(I  +  At  A)  un+1  =  un  (5. 

where  n  denotes  the  time  level,  A  represents  the  spectral  approximation  of  the 
operator  3/3x,  and  un+^  is  obtained  by  inverting  (I  +  At  A)  and  operating 
the  results  on  the  known  quantity  u11.  The  algebraic  equation,  in  a  simple 
form,  can  be  written  as 


Lu  *  f 

(5. 

where 

.  =  (I  +  At  A) 

(5. 

_  n+1 

(5. 

u  =  u 

t-h 

III 

C 

o 

(5. 

Let  L  be  an  approximate  operator  to  L  that  can  be  inverted  easily, 
iterative  scheme  is  suggested  as 


Then  an 


or 


L  (u  -  u  )  ■  -a(Lu  -  f) 
ap  m+1  m  m 


u  -  (I  -  aL_1L)  u  +  aL_1f 
m+1  ap  m  ap 


(5. 


(5.1 


where  a  is  a  parameter  that  can  be  adjusted  to  enhance  the  convergence.  The 
iteration  will  converge  if 

K  =  II  1  -  <L  II 
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Let  B  be  the  central  finite  difference  approximation  of  d/dx.  We  choose 

Ln  =  (I  +  At  B)  (5.12) 

ap 

as  an  approximate  operator  to  L.  The  original  backward  Euler  time-step 
pseudospectral  scheme  can  be  solved  by  the  following  iterative  scheme: 

(I  +  At  B)  u  .  ■  oun  +  (1  -  a)  u  +  At  (B  -  aA)  u  .  (5.13) 

m+1  m  m 


5.2  Convergence  of  Richardson  Iteration 

The  convergence  property  of  the  Richardson  iteration  scheme  presented  above 

is  examined  here.  In  particular,  we  shall  attempt  to  derive  an  expression  for 

the  optimum  parameter  a  in  terms  of  the  spectral  properties  of  the  operators 

L  and  L  . 
ap 


Let  e  =  u  -  u.  Then,  from  Equation  (5.10), 
mm 


e  <  max  ( 1 1-Os | ,  1 1-aS I )  e 

m+i  m 


(5.14) 


where  s  and  S  are  the  bounds  of  the  operator  L  defined  on  a  smooth  test 

ap 

function  x  as 


0  <  a  < 


L_1Lx 

*P 


<  S  < 


N  -*•  »  . 


Optimum  convergence  is  given  by 


opt  s  +  S  ’ 

and  the  corresponding  convergence  rate  is  given  by 

_  S  -  s 


(5.15) 


(5.16) 


Opt  S  ♦  8 


<  1 


(5.17) 


The  bounds  s  and  S  must  be  estimated.  Let  D  =  L  H.  and 

ap 


where 


Then 


DD*Z  -  XZ 
I  I Z 1  |  -  1  . 


.  <  x1?2 

min 


(5.18) 


(5.19) 

(5.20) 


From  Equation  (5.18),  it  is  easy  to  show  that 


T  T 

LL  Z  -  XL  L  Z 
ap  ap 


Hence 


(LZ,  LZ) 

[L  Z  L  Z) 
ap  ap 


Using  the  definitions  given  by  Equations  (5.6)  and  (5.12),  we  have 

_  U  +  At  t(A+AT)Z,  Z)  ♦  At2(AZ,  AZ)} 


{1  ♦  At  [(B+B  )Z,  ZJ  +  At  (BZ,  BZ) } 
For  the  Fourier  method  and  central  difference  scheme, 


T 

A  -  -A 
T 

B  -  -B 


1  +  At  (AZ,  AZ) 
1  +  At2(BZ,  BZ) 


Assume  that  sharp  bounds  exist: 


b2N2  <  (BZ  ,  BZ  )  <  | | B | | 2 
—  max  max  ~ 

(AZ  ,  A Z  )  <  a2N2  <  | |A| |2  . 

max  max  —  — 


1  ♦  (aN  At)2 
_l  ♦  (bN  At)2  _ 


2  I  1/2 


For  the  Fourier  method,  Gottlieb  and  Orszag  (1977)  show  that 


(AZ,  AZ)  _>  (BZ,  BZ)  V"  z  • 


Since  B  is  an  approximation  of  A, 


AZ  .  »  BZ  . 

min  mm 
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and 

then 


s=l  , 


S  -*■  a/b  as  NAt  »  1 


(5.30) 

(5.31) 


Since  Gottlieb  and  Orszag  (1977)  show  that  the  maximum  eigenvalue  of  A  is  it 
times  the  maximum  eigenvalue  of  B,  then  the  scheme  will  not  converge  for 
a  *  1.  From  Equations  (5.16)  and  (5.17),  we  have 


2(b/a) 


opt  1  +  b/a 


1  -  b/a 


opt  1  +  b/i 


(5.32) 

(5.33) 


Note  that  the  optimum  convergence  is  independent  of  N  in  the  limit  NAt  »  1. 
Similar  conclusions  can  be  obtained  for  the  Chebyshev  method. 

The  following  numerical  experiments  have  been  conducted  to  test  the  above 
theory  for  convergence  of  the  iteration  scheme.  The  Chebyshev  collocation 
method  is  used. 

A  smoothed  step  function  is  chosen  as  an  initial  condition.  An  arbitrary 
parameter  a  is  chosen.  Two  computations  are  made  with  different  time  step 
sizes,  At^  and  At^.  From  the  convergence  rate,  one  can  compute  the  bounds  S^ 
and  S2«  By  using  Equation  (5.28),  the  parameters  a  and  b  can  be  computed. 
They  are  then  substituted  into  Equations  (5.28),  (5.16)  and  (5.17)  to  obtain 


the  theoretically  computed  optimum  parameter  a  and  convergence  rate  Kopt< 


These  are  given  for  various  numbers  of  modes  and  different  CFL  numbers  in  the 
right-hand  columns  of  Table  2.  To  verify  the  theoretical  results,  a  is  varied 
for  fixed  N  and  NAt  until  optimal  convergence  is  achieved.  The  experimentally 


obtained  a  _  and  K  are  also  given  in  the  same  table.  The  theoretical  values 
opt  opt  ° 


show  good  agreement  to  the  experimental  values.  The  optimal  convergence  rate 


KQpt  varies  with  the  CFL  number  for  fixed  N  but  is  almost  independent  of  N  for 
NAt  ■  2  as  concluded  above.  However,  slow  convergence  for  high  CFL  numbers  does 


not  imply  higher  computational  costs,  since  fewer  time  steps  are  required  to 
reach  a  specific  time.  The  gain  in  computational  effort  will  be  addressed  later. 
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Table  2.  Numerical  Results  for  Convergence  of  Iteration 

(Chebyschev  Method) 


Numerical  Experiments  I  Theoretical  Values 


1  0.4 

0.80 

0.27 

1  0.86 

1 

1 

| 

0.28 

1  0.8 

0.70 

0.38 

1  0.74 

1 

1 

| 

0.38 

1  1.0 

0.65 

0.43 

1  0.69 

1 

1 

1 

0.43 

1  2.0 

0.40 

0.59 

1  0.51 

1 

1 

0.63 

1  0.4 

0.70 

0.34 

1  0.74 

1 

1 

| 

0.37 

1  0.8 

0.60 

0.43 

1  0.57 

1 

1 

1 

0.43 

1  1.0 

0.50 

0.47 

l  0.51 

1 

1 

1 

0.47 

1  2.0 

0.30 

0.67 

1  0.36 

1 

1 

1 

0.68 

1  0.4 

0.60 

0.37 

1  0.57 

1 

1 

1 

0.38 

1  0.8 

0.50 

0.47 

1  0.41 

1 

1 

1 

0.52 

1  1.0 

0.40 

0.54 

1  0.36 

1 

1 

1 

0.57 

1  2.0 

0.27 

0.69 

1  0.27 

1 

1 

1 

0.69 

r  *•>, 

I 


TR-325/04-85 


-26- 


5.3  Finite  Iteration  Stability 

The  above  section  addressed  the  question  of  convergence  at  each  time 
step.  Here,  we  examine  the  stability  of  the  time-stepping  scheme  with  finite 
Richardson  iterations. 

Consider  again  Equation  (5.5)  and  the  iterative  implicit  scheme. 

Equation  (5.10).  For  a  finite  M-step  iteration,  the  time  stepping  is 


M-l 


u, 


,«  A  ♦  a  y  kV1 

M  o  /  j  ap 


(5.34) 


m=0 


where  uq  is  the  initial  guess  for  the  iteration  and 


K  -  I  -  aL_LL  . 
ap 


(5.35) 


Since 


M-l 


y  k*  -  (i  -  K)_i  (i  -  km) 


(5.36) 


m”0 


and 


L  K®L_1  -  LK^  1  , 

ap  ap 


(5.37) 


Equation  (5.34)  can  be  written  as 


-1  M  -1 

u  “  L  f+K  (u  -L  f)  . 
M  O 


(5.38) 


The  factor  in  the  parentheses  is  simply  the  error  of  the  initial  guess  u 


and  is  denoted  by  e  so  that 

o 


„M 


-  u  +  K"eQ  .  (5.39) 


Let  us  first  choose  an  extreme  example  of  uq  *  0.  Equation  (5.38)  can  then 


be  written  as 


n+1  _  n 

u  *  G  u 


where 


G  *  (I  -  KM)  L_1  . 


(5.40) 

(5.41) 


For  stability,  we  must  have 


llG||  <  1 


(5.42) 
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I 


This  implies  Chat 


I |K| |M  <  p  .  (L)  -  1 
—  mm 


(5.43) 


where  p  .  (L)  is  the  spectral  radius  of  L.  The  estimates  for  the  spectral 
mm 


radius  (Gottlieb  and  Orszag,  1977)  are 


,-1/2 


with 


p  .  (L)  -  !  -  {  CN  .  Chebyshev 

min  '  CN  1  Fourier 


At  -  OOf1). 


(5.44) 


Since,  for  large  NAt,  M K | |  is  independent  of  N  as  discussed  in  the  previous 
section,  Equation  (5.43)  implies  that 

M>0(logN)  (5.45) 

for  convergence. 

Next,  consider  that  the  initial  guess  at  each  time  step  is  given  by  the 
solution  predicted  by  a  forward  Euler  step 


u  ■  (I  -  At  A)  u"  . 
o 

It  can  be  shown  that  the  amplification  factor  is  given  by 

G  -  {KMAt2A2  ♦  1}  L_1  . 

The  stability  criterion  now  gives 

p  .  (L)  -  1 


(5.46) 


(5.47) 


I  1 K |  |M  < 


mm 


At2  p  .  (L)  p  (aV1) 

mm  max 


(5.48) 


But 


p  (A2L_1) 
max 


0(N  )  Chebyshev 
0(N)  Fourier 


(5.49) 


The  stability  requires 


|0(log  N)  Chebyshev 
0(1)  Fourier 


(5.50) 


These  examples  show  that  the  Richardson  implicit  method  can  improve  the 
Chebyshev  method  in  terms  of  computational  efficiency. 

Numerical  experiments  have  been  carried  out  to  substantiate  these  analyses. 
The  Fourier  method  is  applied  to  solve  the  simple  wave  equation  with  the 
periodic  initial  condition  in  (0,  1): 


1  I x-1/2 |  <  1/4 

0  I x-1/2 I  >  1/4 


(5.51) 


The  initial  condition  is  smoothed  by  averaging  the  u  value  of  the  neighboring 
points  twice  using 


u. 

l 


2ui  *  Vi1 


(5.52) 


The  optimal  iteration  parameters  are  chosen  according  to  the  results 

given  in  the  previous  section.  For  each  time  step,  the  iteration  is  terminated 
by  the  condition 

"  <  0  llu  -  u  II  (5.53) 


lu  -  u  I 
m+1  m 


m 


where  ug  is  the  solution  given  by  the  forward  Euler  step  and  o  is  the  tolerance 
parameter.  Table  3  summarizes  the  results  of  the  numerical  experiments  for 
the  Fourier  method.  The  implicit  iterative  scheme  is  stable  for  large  CFL 
numbers  as  compared  to  the  explicit  method.  However,  the  computational  time 
can  be  larger  than  that  using  an  explicit  scheme  due  to  the  large  number  of 
iterations  required  for  each  time  step.  For  the  Fourier  method,  the  iterative 
scheme  does  not  offer  any  advantages  over  an  explicit  scheme. 

We  now  turn  to  the  Chebyshev  method.  Table  4  summarizes  the  numerical 
results.  The  performance  of  the  explicit  scheme  is  very  poor,  as  can  be  seen 
from  the  table.  For  64  modes,  the  method  diverges  even  with  a  CFL  number  of 
0.0125.  The  improvement  in  both  accuracy  and  efficiency  with  an  iterative 
implicit  scheme  is  very  much  evident.  Savings  by  a  factor  of  8  can  be  seen 
for  a  32  mode  calculation.  For  64  modes,  the  iterative  implicit  method  is 
stable  for  a  CFL  number  of  1.0. 
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Table  3.  Fourier  Iterative  Implicit  Method 


T 

1 

1 

N 

I  CFL 

1  L2  Error 

1  Norm 

1  Nell 

CPU/T 

1 

1 

1 

1 

16 

1  0.4 

1 

I  0.54-1 
| 

1 

1 

1 

1 

1 

16 

1  1.0 

1  0.51-1 

1 

1.2 

1 

1 

1 

1 

1 

i 

16* 

1  0.1 

1  0.55-1 

1 

1.3 

1 

1 

1 

1 

1 

32 

1  0.4 

1 

1  0.27-1 

i 

6.1 

1 

1 

1 

1 

1 

i 

32 

1  1.0 

1 

1  0.30-1 

1 

11.0 

1 

1 

1 

1 

1 

i 

32* 

1  0.1 

1 

1  0.28-1 
| 

5.9 

1 

1 

| 

l 

1 

1 

64 

1  0.4 

1 

1  0.13-1 

i 

35.0 

1 

1 

1 

1 

1 

1 

64 

1  1.0 

1  0.19-1 

1 

61.0 

1 

1 

1 

i 

1 

1 

64* 

1  0.1 

1 

1  0.15-1 

1 

26.0 

1 

1 

1 

♦Explicit  Leap-Frog  Scheme 


Table  4.  Chebyshev  Iterative  Implicit  Method 


1 

1 

1 

N 

1 

1  CFL 

|  L2  Error 
!  Norm 

1  Hell 

CPU/T 

i 

i 

1 

1 

16 

1  0.4 

I 

1  0.33-1 

| 

2.1 

i 

i 

1 

1 

1 

16 

1  1.0 

1  0.27-1 

1 

1.6 

i 

i 

i 

1 

1 

1 

16* 

1  0.1 

i 

1  0.33-1 

1 

3.9 

i 

i 

i 

1 

1 

1 

32 

1  0.4 

1 

I  0.17-1 

i 

15.0 

I 

i 

i 

1 

1 

1 

32 

1  1.0 

1 

|  0.14-1 

t 

10.0 

I 

i 

i 

1 

1 

1 

32* 

1  0.05 

i 

1  0.51-1 

l 

31.0 

1 

i 

i 

l 

1 

| 

32* 

1  0.02 

1 

|  0.36-1 

l 

77.0 

i 

i 

i 

1 

1 

1 

64 

1  0.4 

i 

1  0.74-2 

i 

130.0 

i 

i 

i 

1 

1 

1 

64 

1  1.0 

1 

I  0.98-2 

i 

88.0 

1 

i 

i 

1 

1 

1 

64* 

1  0.0125 

1 

|  0.10-7 

1 

400.0 

i 

i 

♦Explicit  Second-Order  Adaras-Bashforth  Scheme 
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6.  Summary  and  Conclusions 

The  present  investigation  attempts  to  construct  a  pseudospectral  scheme 
that  is  highly  accurate  and  competitive  in  computational  efficiency  with 
existing  finite  difference  or  finite  volume  methods.  Instead  of  examining 
simple  model  problems,  we  have  examined  realistic  flows  around  a 
two-dimensional  airfoil.  The  exterior  of  an  airfoil  is  mapped  to  the  interior 
of  a  circle  by  a  conformal  mapping.  We  have  attempted  to  address  several 
questions  in  applying  the  pseudospectral  method  to  transonic  flows  with  shock 
waves.  First,  a  hybrid  scheme  based  on  a  combined  spectral-finite  difference 
method  is  attempted.  Various  means  of  filtering  Gibb's  error  are  examined  for 
their  effects  on  numerical  stability.  Then,  we  attempt  to  construct  a  full 
spectral  scheme  using  spectral  discretization  in  both  spatial  dimensions. 
Finally,  we  have  studied  the  convergence  and  stability  of  an  iterative 
pseudospectral  scheme,  which  could  circumvent  some  of  the  problems  in  the 
construction  of  a  full  spectral  scheme. 

The  following  conclusions  are  made  based  on  the  results  of  our  study: 

(1)  A  hybrid  scheme  using  spectral  decomposition  in  the  direction  along 
the  airfoil  surface  and  a  finite  difference  scheme  in  the  other 
direction  is  capable  of  resolving  shock  waves  in  one  grid. 

(2)  Several  filters  have  been  studied.  A  low-pass  filter  in  the  spectral 
space  is  not  able  to  stabilize  the  computation  without  seriously 
affecting  the  shock  resolution.  An  algebraic  filter  that  averages  the 
flow  variables  around  a  grid  point  is  capable  of  stabilizing  the 
computations  and  maintaining  the  sharpness  of  the  shock  wave. 

(3)  The  filter  being  used  is  stronger  than  the  nonlinear  artificial 
viscosity  expression  given  by  Jameson  et  al.  (1981)  except  near  the 
shock  wave.  Hence,  the  solution  in  the  area  away  from  the  shock  wave 
is  contaminated.  The  drag  count  is  substantially  higher  than  that  in 
a  finite  volume  calculation.  However,  because  of  its  ability  in 
resolving  a  shock  wave,  a  hybrid  scheme  calculation  with  a  coarse  grid 
can  provide  useful  information. 
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(4)  The  time  step  for  an  explicit  hybrid  scheme  is  smaller  than  that  for  a 
finite  volume  calculation,  and  it  takes  a  substantially  larger  number 
of  time  steps  to  achieve  a  steady-state  solution.  An  implicit  scheme 
or  a  multigrid  technique  may  be  required  to  reduce  the  computational 
costs  of  spectral  calculations. 

(5)  The  residue  of  the  scheme  constructed  in  the  present  investigation 
does  not  decrease  with  time.  The  number  of  supersonic  points  in  the 
flow  field  is  taken  as  an  indicator  of  convergence.  Attempts  to  find 
another  form  of  error  norm  have  not  been  successful. 

(6)  An  explicit  full  spectral  scheme  with  a  Chebyshev  polynomial  expansion 
requires  excessive  computing  time  and  is  not  competitive  with  a  finite 
volume  calculation  using  a  dense  grid.  A  polynomial  subtraction 
Fourier  method  also  does  not  show  much  promise  as  a  superior  method  to 
the  finite  volume  method. 

(7)  A  Richardson  iterative  spectral  implicit  scheme  has  been  analyzed  for 
a  simple  wave  equation.  The  analyses  and  the  accompanying  numerical 
experiments  show  that  there  exits  an  optimum  convergence  rate  for  the 
iteration.  In  order  for  the  computations  to  be  stable  for  a  finite 
iteration  scheme,  a  certain  minimum  number  of  iterations  is  required 
for  each  time  step.  As  a  result,  the  Fourier  implicit  iterative 
scheme  does  not  offer  any  savings  in  computational  cost  for  a  steady- 
state  problem.  Using  a  Chebyshev  method,  savings  by  a  factor  of  N  in 
operation  count  can  be  achieved  by  using  an  iterative  implicit  method 
rather  than  an  explicit  one. 
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Figure  4.  Finite  Volume  Calculation — Distribution  of  Pressure  Coefficient 
on  64x16  Grid  (Mo*  0.63;  a  *  2  degrees;  CL  -  0.3419; 

CD  -  0.0018;  CM  -  -0.0078) 
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Figure  6.  Hybrid  Scheme — Initial  Diatribution  of  Presaure  Coefficient 
on  64x16  Grid  (M*  ■  0.7;  a  ■  2  degrees;  *  0.022; 

CD  ■  0.002) 
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Figure  7.  Hybrid  Scheme — DistribuCion  of  Pressure  Coefficient  on  64x16 
Grid  After  800  Time  Steps  (M<*>  *  0.7;  a  *  2  degrees; 

CL  -  0.411;  CD  -  0.019) 
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Figure  8.  Hybrid  Scheme — Distribution  of  Pressure  Coefficient  on  64x16 
Grid  After  2400  Time  Steps  (Moo  ■  0.7;  a  •  2  degrees; 

CL  -  0.323;  Cj)  -  0.021) 
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Figure  10.  Finice  Volume  Calculation — Distribution  of  Pressure  Coefficient 
on  64x16  Grid  (M®  ■  0.7;  a  ■  2  degrees;  Cl  ■  0.3699; 

CD  -  0.0131;  CM  -  -0.0118) 
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Figure  13.  Chebyshev  Method — Pressure  Waves  for  One-Dimensional  Problem 
Using  33  Grid  Points 
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